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We develop our recently proposed lattice-Boltzmann method for the non-equilibrium dynamics of 
amphiphilic fluids Jl] . Our method maintains an orientational degree of freedom for the amphiphilic 
species and models fluid interactions at a microscopic level by introducing self-consistent mean-field 
forces between the particles into the lattice-Boltzmann dynamics, in a way that is consistent with 
kinetic theory. We present the results of extensive simulations in two dimensions which demonstrate 
the ability of our model to capture the correct phenomenology of binary and ternary amphiphilic 
fluids. 



I. INTRODUCTION 



The lattice-Boltzmann (LB) method is a mesoscopic approach to the study of the dynamics of fluids, intermediate 
between macroscopic models, based on Navier-Stokes equations, and microscopic approaches, based on molecular 
dynamics (MD) simulations. In this method, time, velocity and space are all discretized and a single-particle distri- 
bution function is utilized to describe the time evolution of an ensemble of molecules having a discrete set of possible 
velocities ^,[| . The computationally demanding tracking of individual molecules is thus avoided and at the same time 
the ensemble-averaged distribution function retains much of the microscopic information. Macroscopic or hydrody- 
namic effects naturally emerge from mesoscale lattice-Boltzmann dynamics, provided that the LB equation possesses 
the correct and necessary conservation laws and symmetries [0. Historically, the lattice-Boltzmann model was first 
derived from its predecessor-the lattice-gas automata (LGA) but it has been shown recently that the model 

can also be derived from kinetic theory by discretization of the Boltzmann equation in velocity, space and time [^,^| . 

An important and promising area of application of the lattice-Boltzmann method is for modeling the dynamics of 
multicomponent fluids. There is in particular much fundamental and technological interest in modeling amphiphilic 
systems; i.e. fluids comprising one or two immiscible phases (such as oil and water) and an amphiphilic species 
(surfactant). The addition of surfactant to a binary immiscible fluid can result in complex structures on a mesoscopic 
length scale, including lamellar, micelles and microemulsion phases The non-equilibrium dynamics and hy- 

drodynamics of these systems are difficult to simulate using a continuum approach based on Navier-Stokes equations. 
One major difficulty these methods face is the existence of complicated interfaces between fluid components which can 
undergo topological changes; another is that it is far from clear how to formulate a hydrodynamic description of these 
fluids. Microscopic approaches, based on molecular dynamics, are able to deal with amphiphilic fluids, but they are 
still computationally too demanding to investigate large-time dynamics (times scales accessible to MD are typically of 
the order of nanoseconds [[[o| whereas important time scales in such complex fluids vary between 10~ n — 10 4 seconds) 
and the extended spatial structures involved in many problems of interest. 

The formation of fluid interfaces is a consequence of interaction between the molecules of the fluids []l3| , and the LB 
models for multicomponent fluids come in two varieties, depending on the way interactions between fluid components 
are incorporated. In the free-energy approach [13-16] the starting point is an Ansatz for the free-energy functional of 
the complex interacting fluid and the approach to equilibrium is assumed to be governed by a Ginzburg-Landau or 
Cahn-Hillard equation or some suitable generalization thereof. A LB collision operator is then constructed that gives 
rise to the desired evolution equation in the hydrodynamic limit. The most attractive feature of this model is that the 
equilibrium state is fixed a priori by the choice of free-energy functional. This guarantees thermodynamic consistency 
at equilibrium but limits the description of non-equilibrium dynamics. For example, the one-component free energy 
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model is not Galilean invariant |l7fl . Furthermore, the method relies on the ability to postulate a suitable form for the 
free energy, which might not be always available. An alternative approach, which is closer in spirit to kinetic theory, 
is to add interactions between fluid components by introducing intermolecular forces between fluid particles pi| , p5[ . 
Very recently we proposed such a lattice-Boltzmann model for amphiphilic fluids JjJ which is based on introducing the 
dynamics and interactions of amphiphilic molecules into the multicomponent lattice-Boltzmann scheme of Shan and 
Chen In the resulting ternary vector model, which was inspired by the lattice-gas automata (LGA) formalism 
of Boghosian Coveney and Emerton 26 1 , amphiphilic species possess both translational and orientational degrees of 
freedom, a feature which is crucial for capturing much of the above-described phenomenology of amphiphilic systems. 

The principal aim of the present paper is to further develop and explore this model and to assess and establish 
its general validity for modeling the dynamics of amphiphilic systems. In Section II we give a detailed description 
of the model and show that with slight modification, the mean-field treatment of intermolecular interactions in the 
model |i|,[l8| can be made consistent with the mean-field kinetic theory of interacting fluids, based on a set of coupled 
Boltzmann equations. In section III we use this modified model to simulate the non-equilibrium dynamics of binary 
immiscible and binary and ternary amphiphilic fluids in two-dimensions. These simulations demonstrate that the 
model exhibits the correct dynamical behavior in the case of critical quenches of binary fluids and it is capable of 
describing in a consistent way the phenomenology of various experimentally observed self-assembling structures, such 
as lamellae, droplet and sponge microemulsion phases, the arrest of separation of immiscible oil-water phases when 
enough surfactant is present in the system, and the formation of the lamellar phases in binary water-surfactant systems. 
Finally, in section V we draw some conclusions from this work and discuss the prospects of future developments of 
the model. 



II. LATTICE-BOLTZMANN MODEL FOR INTERACTING FLUID MIXTURES 



The lattice-Boltzmann scheme for a fluid mixture having S components is defined in terms of the distribution 
functions f?(x,t) belonging to the particles of component er of the fluid mixture (e.g. a = water or oil in a binary 
system). 

In the BGK approximation |^0| of the collision operator these equations are |^|| 

/f (x + Cl At,t + At) - f? (x,f) = - fr ~. fi " At (1) 

-v 

where Ci (i = 0, 1 . . . M) are a set of discrete velocity vectors on a regular lattice, x is a point on the underlying 
spatial grid, f^ eq ^ is a local equilibrium distribution function, A CT is the relaxation time for component a and At is 
the timestep. 

The number density n a and the mean velocity u CT of component a are given by 

n CT (x,i) = £/f, (2) 



(3) 



with p a = m r, n a the density and m a the mass of species a. 

The equilibrium distribution function should be chosen in such a way that the BGK collision operator f2° 

~(fi ~ fi^ B ^)/^a- locally conserves the mass of each species a 



(4) 



We also require that the momentum of all components together should be conserved locally 

^m^CiQ* = 0. 



(5) 



A common choice which satisfies the above constraints is l2llJ2z 



rcr(eq) 



2cl 



(6) 
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where u>i is a weight factor and c s is the speed of sound, both of which are determined entirely by the choice of lattice. 
For example, in the so-called D3Q19 lattice j2^] (19 velocity vectors in three dimensions) c s = and u>i = 1/3, 

1/18 and 1/36 for ei = 0, = 1 and e, = 2, respectively. 

The requirement that the total velocity of all components together should be conserved at each lattice point results 
in the following form for the common velocity u 

n = ^r^- (7) 

In the limit of small Mach numbers, i.e. u/c s -C 1, where (6) is positive, the above choice of /f^ 9 ' results in 
Navier-Stokes hydrodynamics, where each fluid component a has kinematic viscosity v a = (X a /At — l/2)c? s |2l[| . 



Mean-field treatment of interactions 



In order to describe immiscibility effects, some form of repulsive interaction between fluid components must be 
introduced into the LB equations. In the Shan-Chen scheme ^M , coupling between fluid components is achieved by 
including intcrmolccular interactions in the model. These interactions are modeled as a self-consistently generated 
mean-field force. Allowing only homogeneous isotropic interactions between nearest-neighbors, the mean-field force 
F" acting on particles of component a is given by 



F ff (x, t) = -VAx, t) J2 E S-^( x '> *)(*' - x) 



(8) 



where ^ a = ip' y (p a (x, t)) is the so-called effective mass, which can have a general form for modeling various types of 
fluids (in the present study we have taken tp — p a ) and g a „ is a force coupling constant, whose magnitude controls 
the strength of the interaction between component a and a and whose sign determines whether this interaction is 
repulsive {g aS > 0) or attractive (g aS < 0). Shan and Chen incorporated the above force term in the LB dynamics 
by adding an increment ffq] 



Su 17 = — r a At 
p° 



(9) 



to the velocity u which enters the equilibrium distribution function (6) in the BGK collision operator. 

With the choice ip a = p a the above treatment of intermolecular interactions is consistent with kinetic theories of 
interacting fluids, in which non-local interactions are described by including a self-consistently generated mean-field 
force term in the Boltzmann equation while local interactions are treated as true collisions. Indeed, kinetic equations 
of this type have been successfully applied to describe gas-gas separation into two phases and, more recently, 
spinodal decomposition in binary fluids p7| . However, the way the mean-field force enters the Shan-Chen model is 
not entirely consistent with the Boltzmann equation in the presence of a force [^8[^9) . The explicit expression for the 
force term in the Shan-Chen model can be obtained by substituting u + i5u into the equilibrium distribution function. 
This yields 



/?(x + CiAt,t + At)-/?(x,t) 



tf-fl 



eg) 



- UJiTl 



1 



(C 



(Cj-U) 

■a CT ) 2 



.a a At 
■ At 2 . 



(10) 



where a CT = F cr / ( o <T is force per unit mass and we have introduced the dimensionless relaxation time t„ = \ a /At. 
As was shown by Luo J28| and by Martys et al. |^9| , by omitting the second order term in & a the above equation 
can be made consistent with a systematic derivation of the lattice-Boltzmann model in the presence of the force a CT , 
as obtained from discretization of the corresponding Boltzmann equation, in the low Mach-number regime. In our 
scheme we omit this non-linear term when calculating the effect of the self-consistently generated mean-field force, so 
that the resulting lattice-Boltzmann scheme is consistent with the underlying Boltzmann equation in the presence of 
a force. Our final lattice-Boltzmann equation is then given by 



/f(x + Ci At,t + At)-/f(x,t) 



n-fi 



(Cj.u) 



SL„At 



(11) 
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In the case of a single-component interacting fluid rewriting the force term in terms of distribution functions nf 
leads to the following lattice-Boltzmann equation 



/f(x + c i At,t + At)-/f(x ) i) = - 



f[ f, 



<Ti(eq) 



(12) 



where A^ 7 is given by 



A CTCT 



r(ct ~Cj) 



.a a At 



(13) 



Thus the net effect of the force term is to introduce off-diagonal matrix elements in the BGK collision operators. 
These matrix elements describe transitions from state Cj to state Cj, and vice versa, due to the drift of particles 
under the influence of the force during the time At. In the case of a fluid mixture the BGK collision operator couples 
different fluid components even in the absence of the mean-filed force, through the common velocity u which enters 



cr(eg) 



The mean-field force introduces an additional non-local coupling between fluid particles and in this case Eq. 



(11) can be rewritten as 



f?(x + Ci At,t + At)-f[(x,t) 



ft - f 



(14) 



where A^ CT are matrix elements of the cross-collision operator 



2 {.^aa^i ^ererC?) 



C..C; 



.a ff Ai 



(15) 



with 



rr 



L-^O TV, 



(16) 



The above treatment of intermolecular interactions does not take into account strong repulsive interactions between 
molecules at short distance which prevent them to overlay (the excluded volume effect [Q ) . This effect is important 
in calculating transport coefficients of dense gases and fluids but we do not expect it to significantly alter the 
dynamics of microemulsion, which is the main focus of our present work. The excluded volume effect can be introduced 
in the model by adding Enskog corrections to the BGK collision operator J2q] . 



A. Modeling surfactant dynamics and interactions 



A surfactant (or amphiphilic) molecule usually possesses two different fragments, each having an affinity for one 
of the two immiscible components. In the case of a prototype mixture of immiscible fluids like oil and water, a 
typical surfactant molecule would be an amphiphile that has a hydrophilic head preferring to be in contact with water 
molecules, and a hydrophobic tail preferring to be surrounded by oil. Under equilibrium conditions the surfactants 
are predominantly adsorbed at the oil-water interface, hence effectively screening the two-body repulsive interaction 
between immiscible species. In our scheme these essential characteristics of amphiphilic molecules are modeled 
by introducing a separate amphiphilic species into the Shan-Chen model which possesses both translational and 
rotational degrees of freedom. This is achieved by modeling these species as idealized point dipoles. Consequently, 
the interaction of the amphiphilic molecules with each other and with non-amphiphilic molecules (like oil and water) 
depends not only on their relative distance but also on their dipolar orientation p6fl . Furthermore, a full description 
of the dynamics of amphiphilic molecules requires a description of the propagation of each single-particle distribution 
function as well as the time evolution of their dipole vectors. Making the physically reasonable assumption that the 
dipole moment carried by amphiphilic molecules is independent of their velocity, an average dipole vector d(x, t) 
is introduced at each site x representing the orientation of any amphiphile present there. The propagation of the 
amphiphilic molecules is then described in our scheme by a set of BGK-likc equations, one for the distribution function 
// and one for the relaxation of the average dipole vector d(x, t) to its local equilibrium orientation [0: 
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//(x + CiAt,t + A*)-//(x,t) 



s(eg) 



(17) 



d(x, t + At) = d(x, t) [d(x, t) - d eq (x, t)] . 

Td 



(18) 



In the above equations /? (x, t) is the distribution function of surfactant molecules and r s and Td are dimensionless 
relaxation times representing the rates of relaxation of //(x, i) and d(x, t) to their local equilibrium distributions 

f^ eq ^ and d e9 (x, t), respectively. The collision operators Af } s and AfJ 5 have the same form as (15) and describe the 
effect of the total mean-field force experienced by surfactant particles due to their interactions with non-amphiphilic 
particles and with other amphiphilic particles. We shall specify these forces shortly. Finally, d(x, t) is obtained from 
the following equation of motion 



5 (x, t)d(x, t) = /*( x - c 4 At)d(x - CiAt, t). 



(19) 



The equilibrium distribution function f^ eq ^ is chosen to have the same form as (||) and, in analogy with the Weiss 
molecular field theory of magnetism p0[ , the equilibrium distribution d eq (x, t) is obtained self-consistently from the 
Boltzmann distribution as 



d e9 (x,f) = d 



(20) 



where dil is an element of solid angle, v is a unit vector and U is the potential energy of a dipole in the mean field 
generated by dipolar amphiphiles and non-amphiphilic molecules: 



U ■■ 



-v • (h s 



(21) 



In ( pi]) h s and h c are the mean fields resulting from dipole-dipole interactions and the interaction between surfactant 
dipoles and non-amphiphilic molecules, respectively. In Eq. (^0|) both dg and the inverse temperature /3, are 
independent parameters. Performing the integral yields for the equilibrium distribution 



d 



coth (J3h) 



0h 



in 3D, and 



d e? = d 



(22) 



(23) 



in 2D, where h = h/h and Iq and I\ are the zero and the first order modified Bessel functions [Q, respectively. 
The mean-field generated by water and oil molecules is given by 



h c (x, t) = ^ E ™ CT ( X + c * 



(24) 



where e„ is the "charge" for various molecular components (which may take its values from the set {—1,0, 1}). In 
the present simulations we take e = 1 for water molecules and e — — 1 for oil molecules. Similarly, the mean-field 
generated by other surfactant molecules is given by 



h s (x,£) = ^[^nf(x + CjAt,i)0j • d*(x + CjAi, t) + n|(x, i)d,(x, t)}, 

i j^O 



(25) 



where 



is the traceless second-rank tensor, and D is dimension of the lattice. 
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Finally, we specify the form of interaction between amphiphilic molecules with water and oil particles and between 
amphiphilic molecules themselves. These are obtained from Eq. (|^) by treating each amphiphilic molecule as a pair 
of water and oil molecules displaced by a distance d(x, i) from each other and performing a Taylor expansion in d 
in the resulting expression for the total force jD. Assuming that the dipole head and tail have equal and opposite 
charges e = ±1 and only nearest-neighbor interactions considered, the additional forces are given by 



F ff < s (x, t) = -2^ CT (x, t)g as d(x + c, At) • (I - -^D)V> s (x + c.At, t) (26) 

ifto c * 



F^ c (x, t) = 2^ s (x, f)d(x, t) ■ £(i - yB)f(x + a, t) (27) 



and 



AT) c c 

F^(x,t) = 2-^ s (x)^{d(x + Ci At,t)d(x,i) : [I--^D] Ci 

c i c * 

+ [d(x + c t At, t)d(x, t) + d(x, t)d(x + c,At, t)] ■ c ( }^ s (x + c; + At, t) (28) 

In the above equations F"^ is the force acting on non-amphiphilic particles a (water and oil) due to amphiphile 
dipolcs, F s,c is the force acting on amphiphilic molecules due to all non-amphiphilic particles and F s,s is the force 
among amphiphilic molecules themselves. The coupling constants g GS and g ss determine, respectively, the strength of 
interaction between water/oil molecules and surfactant molecules, and among surfactant molecules themselves. The 
coupling coefficient g ss should be chosen negative if we wish to model attraction between two amphiphile heads or 
tails, and repulsion between a head and a tail. 



III. SIMULATIONS 



As mentioned earlier the principal aim of the present work is to assess the ability of our model to reproduce the basic 
properties of self-assembling amphiphilic fluids. For this reason we choose not to explore the entire parameter space 
of the model in the present study but rather find through a limited search in this space a canonical set of parameters 
which allows us to describe generic behavior in a consistent way. Another important consideration in choosing the 
parameters relates to the numerical instabilities which, for a given concentration of oil, water and surfactant, were 
found to occur when the force coupling constant g aS , g as and g ss , or the mean densities, were increased beyond certain 
threshold values. We found that these instabilities occur when the forcing terms cause the right-hand side of Eqs. (14) 
and (17) to become negative and are caused by a combination of the mean-field treatment of interparticle interactions 
and the restriction of the lattice-Boltzmann scheme to low Mach numbers. With these considerations in mind, and 
after a restricted search in the parameter space of the model we arrived at the following set of canonical parameters 
which, unless stated otherwise, are used throughout our simulations (the timestep At is set to 1 throughout) g aa = 0, 
g aS = 0.03, g as — —0.01, g ss = 0.01, r a = t s = 1, = 2, m a = 1, m s = 2, do = 1 and (3 = 10. All calculations were 
performed on a 256 2 lattice subject to periodic boundary conditions. In the case of the binary oil-water systems we 
also performed additional calculations on 512 2 and 128 2 lattices in order to check for finite size effects. 

The CPU time and memory requirements of our algorithm both scale ~ L D , where L is the linear dimension 
of the lattice. In two dimensions we were able to study adequate system sizes using typical workstations. For 
example, simulation of a ternary system on a 256 2 lattice on a Silicon Graphics 250 MHz processor workstation takes 
just under 8 hours to reach 3000 timesteps (switching off the subroutines which perform the surfactant dynamics 
calculations reduces the CPU time by 50%.) For the present calculations we used a serial implementation of our 
lattice-Boltzmann algorithm. In three dimensions, however, the serial algorithm quickly becomes prohibitive in terms 
of computer memory for moderately sized systems. Fortunately, an important feature of LB is its intrinsically parallel 
structure and we have implemented a parallel version of our algorithm [fl9[ which allows us to perform very large- 
scale ZD simulations on massively parallel platforms. We plan to report on the three dimensional model in future 
publications. 



A. Binary Oil- Water System 



The dynamics of phase separation in a binary mixture, following a thermal quench into the unstable coexistence 
regime, proceeds by spinodal decomposition. Immediately after the quench, small domains, with local concentrations 



G 



roughly corresponding to that of the two pure immiscible phases, spontaneously form and grow and finally result in 
complete phase separation [|l2| . To simulate phase separation we set up a simulation with equal quantities of water 
and oil with average densities 0.5, and no surfactant present (a "critical quench"). This choice of average densities 
ensures that we are well within the immiscibility region of the model. The initial condition of these simulations is 
a uniform mixture of the two fluids with small random fluctuations in the uniform densities. These fluctuations are 
necessary to move the system out of a metastable uniform state, in which the mean-field forces are identically zero. 
The force term is initially set to zero and ff are set to f^ eq ' calculated from if and u = 0, (Eq. (7)). 

As can be seen from Fig. 1, immediately after the quench small domains spontaneously start to form. As time 
evolves, sharp interfaces develop between the regions associated with each phase, and the branchlike structures which 
were formed at the earlier stage of simulations coarsen. The growth of domains continues and, if left to run for a 
large enough time, the system would eventually reach the completely separated state of two distinct layers of fluid. 
Phase-separation experiments typically measure the structure factor, S*(k, t), which contains information on the time 
evolution of the various length scales present in the system. It is defined as the Fourier transform of the density-density 
correlation function. For the discrete systems we are studying, we consider equivalently 



5>(x,<)-«^)K 



(29) 



where k = (2w/L)(mi + nj), m,n — 1,2,... L; here L is the linear lattice size, N — L 2 is the total number of grid 
points in the system, q(x.,i) = n water (x,t) — n (x,t) is the total color order parameter at grid point x and time t 
and q(t) is the spatial average of q(x,t) at time t. S(\c, t) is further smoothed by averaging over an entire shell in k 
space to obtain the circularly-averaged structure factor 

fif(M) = £*?M (30) 

where the sum is over a circular shell defined by (n — \) < < (n + |) and the cutoff wavevector k c has 
the maximum possible value which is compatible with the periodicity in k space. We use the first moment of the 
circularly-averaged structure factor as a measure of the characteristic length scale of the system 

m ~ E k s(k,t)- (31) 

The characteristic size of the domains is then given by R(t) = 2ir/k(t). 

Fig. 2 displays the temporal evolution of the circularly-averaged structure factor obtained for a critical quench in 
a 512 2 system. At early times in the simulations (t < 4000 timesteps) we observed that the amplitude of the peak 
in the structure factor increases without the position of the peak changing in time. This behavior is indicative of 
initial sharpening of the domains, as the amplitude of the peak in S(k, t) is proportional to the domain mass within a 
characteristic domain size. As time evolves further, the peak of S(k,t) shifts to smaller wavenumbers and its height 
increases, indicating the coarsening of the domains. At intermediate times we also observe the appearance of a second 
peak in S(k, t), which later on becomes a shoulder of the main peak and then merges with it. After an initial transient 
regime, the dynamics of phase separation is often characterized by a single time scale. This length scale is usually 
described by a power law behavior R(t) = t a . Simple dimensional analysis of the hydrodynamical evolution equations 
in 2D yields for the domain growth exponent a = 2/3, when the binary system is in the inertial hydrodynamic regime 

us 

In Fig. 3 the time evolution of R(t) is shown as obtained from our lattice-Boltzmann simulations with system size 
512 2 . Finite-size effects in this quantity are known to become important when R(t) is comparable to the linear size 
of the lattice L. We checked these effects by performing additional calculations of R(t) for systems size 128 2 and 
256 2 . By comparing the results we deduced that finite-size errors in our 2D simulations become important when 
R(t) > L/5. The 3D lattice-Boltzmann simulations of Kendon et al. (34|] seem to indicate a somewhat larger value of 
R(t) beyond which these errors become pronounced. This might indicate that finite-size effects are more significant in 
2D. Discarding both the early-time transient regime and the late-time regime where finite-size effects are pronounced, 
we found that the late-time behavior of R(t) in our simulations is R(t) ~ ^o.66±o.oi^ m g 00c j agreement with the above- 
mentioned scaling arguments and previous lattice-gas |36f and lattice-Boltzmann [fh0[ simulations of phase separation 
in 2D. We note, however, that this result is only a first qualitative study of the dynamics of phase separation within 
our model and more work is needed in order to unambiguously identify different scaling regimes within the parameter 
space of the model. 
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B. Microemulsion phases: oil droplets in water and sponge phase 



We used our model to simulate the different ternary microemulsion phases that are possible in 2D, namely the oil- 
in-water droplet and sponge phases. In experimental systems the two distinct microemulsion phases will form when 
the appropriate concentrations of oil, water and surfactant are present in the system below the critical temperature. 
The oil-in-water droplet phase typically consists of finely divided spherical regions of oil, with stabilizing monolayers of 
surfactant surrounding them, embedded within a continuous water background. If one increases the relative amount of 
oil in the system and there is sufficient amphiphile present, one observes the formation of mutually percolating tubular 
regions of oil in water, with a monolayer of surfactant sitting at the interface. In both these cases, the equilibrium 
state does not correspond to complete separation of immiscible oil and water regions, but rather to complex structures 
with very different characteristic length scales that form as the result of the presence of amphiphile jjjj. 




In order to reproduce the oil-in-water droplet phase, we set up a simulation with a uniform mixture of oil, water and 
surfactant with small density fluctuations as our initial configuration. The average densities of oil, water and surfactant 
are 0.42,0.6 and 0.1 respectively. Fig. 4 displays the results. We see the rapid formation of many small oil-in-water 
droplets, whose size initially increases slightly, but not without limit. This is characteristic of the experimental droplet 
phase. It occurs because the free energy penalty associated with the existence of many oil-water interfaces in this 
phase, as compared to the complete oil-water phase separation, is compensated by the gain in free energy due to 
adsorption of surfactant dipoles at these interfaces. If coarsening of oil droplets was to continue the interface area for 
adsorption of surfactant dipoles would decrease resulting in an increase in the amount of surfactant in bulk water or 
oil. The free energy penalty for bulk surfactant prevents this to happen. The concentration of the surfactant is not 
visible in Fig. 4, but is high at the interface and low in oil-rich and water-rich regions, with the surfactant dipoles at 
the interface pointing on average from water-rich region towards the oil-reach droplets (see also Fig. 11 and Fig. 2 in 



In order to further quantify the result we show in Fig. 5 the time evolution of the circularly-averaged structure 
factor S(k, t). Once again we observe the formation of a distinct peak in the structure factor, indicating the sharpening 
of the domains. Interestingly, this happens much faster than in the case of the binary system. The presence of the 
surfactant seems to accelerate domain formation in the early stage of phase separation, an effect which has also been 
seen in the hybrid Ginzburg-Landau simulations of Kawakatsu et al. |35| , and should be detectable experimentally. 
As time evolves the peak in S(k, t) increases in height and shifts further towards smaller values of the wavelength, 
indicating the growth of droplets. From at least timestep 5000 onwards there appears to be a negligible amount of 
further movement of the position of the peak, indicating that droplets have reached a maximum size and will grow 
no more. We observed small oscillations in the intensity of the peak in S(k, t) which suggest that the characteristic 
domain mass fluctuates in time. 

To investigate the ability of the model to access the sponge phase, we set the average densities of water and oil 
both equal to 0.5 while keeping the average density of surfactant at the same value as before. The results are shown 
in Fig. 6. Starting once again from a perturbed uniform mixture of fluids, this time we observe the growth of an 
interconnected network of tubular regions of oil and water. The width of the oil and water regions grows in size 
up to about 4000 timesteps. During this time the surfactant particles, which were initially distributed uniformly in 
the system, concentrate around the various oil-water interfaces. Beyond this stage the system changes very little, 
indicating that the observed sponge phase is stable. This is also confirmed by our result for the circularly-averaged 
structure factor shown in Fig. 7. 



To further investigate the effect of surfactant on domain growth dynamics of the sponge-phase we performed 
additional simulations in which we kept the average density of oil and water fixed at 0.5 while gradually increas- 
ing the amount of surfactant in the system. The average densities of surfactant used in these simulations are 
0.05,0.1,0.15,0.20,0.30. We analyze the effect of varying the surfactant concentration using the domain size R{t) 
as a quantitative measure. The domain size is calculated from the circularly-averaged structure factor, as described 
in Section II. The results are summarized in Fig. 8 where the domain sizes are plotted against time, and as a function 
of increasing surfactant concentration. The presence of the surfactant significantly retards the growth of the domains 
and it can be clearly seen that for surfactant concentrations larger than 0.05, the domain size reaches saturation. 
Following Boghosian, Coveney and Emerton |36] ], we analyzed the time-dependence of domain growth in terms of a 
stretched exponential form 




0)- 



C. The effect of surfactant on domain growth dynamics 



R{t) = R, 



'OO 



aexp(—ct d ) 



(32) 
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where i?oo, a, c and d are adjustable parameters, which are determined by a least-square fit to our data. As can be seen 
from Fig. 8 this form fits our results extremely well across the full time scale of the simulations and for all surfactant 
concentrations considered. We also investigated a fit of the logarithmic form R(t) = a + 6(lni) c , which describes the 
phase-separation of binary alloys in the presence of impurities [ j37| . Obviously, this form is unable to describe the 
late time saturation of the domain size and we found that the root-mean-square errors using this form to describe 
the early times of domain formation are also an order of magnitude larger than that of the stretched exponential 
form. In the case with average surfactant density equal to 0.05 the domain size does not saturate, nevertheless Eq. 
(32) provides a better fit to the slow growth of the domain size than the logarithmic form for the time interval we 
considered (0 < t < 12000). LGA simulations ]36|,[58| indicate, however, that the logarithmic from is a good fit for 
describing the dynamics of "metastable" (i.e. long-lived) emulsions which do eventually phase separate. 

At late times in these simulations, small but persistent oscillations in R(t) can be seen in Fig. 8, which are absent 
in the case of binary systems. These oscillations have also been observed in LGA simulations, both in 2D and 
in 3-D p8| . Their presence in our lattice-Boltzmann simulations confirms that these oscillations are caused by the 
additional dynamics which the presence of amphiphile introduces into the system: they are not a consequence of 
statistical fluctuations inherent in LGA. 

Finally, we investigate the relationship between the final size of the domains and the surfactant concentration. 
Laradji et al. fflfl used simplified MD simulations to study the dynamics of ternary oil, water and surfactant system 
and found the final domain size Rf is inversely proportional to the average concentration of surfactant molecules n s : 

R f ~ — . (33) 

n s 

In the deep quenches with no system fluctuations performed by Laradji et al surfactant molecules entirely reside at 
the interfaces. This was not the case in the lattice-gas simulations of Emerton et al. |36|] in which a certain amount of 
surfactant existed as monomer in bulk oil or water. After subtracting away from n s a background monomer density, 
these authors also found a linear relationship between Rf and l/n s . A similar situation to lattice-gas simulations 
occurs in our simulations where a significant fraction of surfactant exists as monomers in bulk oil or water. However, 
we found that even without correcting for the background monomer density of surfactant Eq. ( |3^) gives a very good 
description of the relationship between Rf and n s . This is shown in Fig. 9 where Rf is plotted versus the inverse of 
n s . The condition that all surfactant molecules should be at the interface does not therefore seem necessary for ( |33| ) 
to hold, as long as the surfactant molecules are mainly concentrated at the interface. 



D. Ternary phase: lamellae 



Next we use our model to investigate the stability of a lamellar structure, which is composed of alternating layers of 
oil-rich and water-rich phases separated by a layer of surfactant. A preliminary discussion was given in jjj. We look at 
the system with and without surfactant present in order for a critical comparison to be made. A similar investigation 
of the lamellar structure in 2D and in 3D has been performed previously using LGA [ p6|j39| . In a similar way to these 
simulations, we set up the initial configuration of the system as layers of pure oil and pure water eight sites wide such 
that each species has an average density 1.0. 

Snapshots of our simulations for the binary case are shown in Fig. 10 for timesteps t < 9000. As can be seen 
from this figure, the initial layered structure becomes less sharply defined as time evolves but the lamellar structure 
remains intact and does not evolve to complete phase separation. By letting the simulations run for much longer 
times (30000 timesteps), we checked that the lamellar structure is indeed the final equilibrium state of the system and 
we are not observing a metastable state with long equilibration time. Wc also examined the stability of the structure 
against small fluctuations in the density and found the lamellar structure to remain stable. These results are in sharp 
contrast with the previous LGA simulations in which, starting from a layered structure, complete phase separation 
was observed |26|] , As pointed out in |lj, in an infinite two-dimensional system with finite-temperature fluctuations, 
one expects lamellar structures to be unstable, due to Peierls instabilities Q. The Peierls theorem, however, does 
not make any statement about the stability of such structures in finite 2D systems. The stability of the lamellar 
structure seen in our lattice-Boltzmann simulations and its instability in LGA simulations thus provides numerical 
evidence that the Peierls mechanism is also capable of destabilizing periodic layered structures in finite systems. 

Next we examine the effect of surfactant on the layered structure by setting up a simulation where there is, in 
addition to water and oil layers, a single layer of surfactant at each of the oil-water interfaces. The average density 
of surfactant in each monolayer was 1.0 and the initial condition for the surfactant dipole vectors at each site was 
d(x, 0) = 0. We found that the final state of the system is, once again, lamellar. The presence of surfactant, however, 
causes the water-oil interfaces to remain sharper than in the previous simulations. This effect is shown in Fig. 11 



9 



(top panel) where the initial and final state of the color order parameter (averaged over the y direction, the direction 
perpendicular to layers) are plotted along the rc-axis for both sets of simulations. 

Due to the symmetry of the system d y , the component of the dipole vector in the y— direction, does not evolve in 
time from the initial condition d y = 0. Under, the influence of the field set up by the color gradient, however, d x 
does evolve in time. A plot of the final state of d x , averaged in the y direction, is shown in figure 11 (lower panel). 
It can be seen that the surfactant directors have their largest values around the interfaces, i.e. at the points where 
the color order parameter changes sign. It is interesting, albeit expected, that the surfactant dipoles alter alignment, 
always pointing from water-rich to oil-rich layers. The reason for this behavior is that, neglecting surfactant-surfactant 
interactions, the direction of the equilibrium dipole vectors in our model is determined by the gradient of the color 
order parameter, as can be seen by expanding Eq. (24) in a Taylor series in the ratio of |cj to the color gradient scale 
length. For the lamellar structure, the color order parameter changes only in the x— directions with its slope changing 
sign alternately, in passing from a water- rich layer to an oil-rich layer, as can be seen in Fig. 11 the surfactant dipole 
vectors flip direction every time this happens. 

E. Self-assembly in mixtures of water and surfactant 

Repeating the simulations performed for the ternary mixtures but setting the average density of oil (or water) to zero 
gives results for binary water and surfactant system. We kept the average density of water fixed at 0.5 and performed 
two simulations, one with a high surfactant density n s = 0.25 and one with a low surfactant density n s = 0.071. 
Snapshots of the simulation for the case of high surfactant concentration are shown in Fig. 12. It can be seen that 
starting from a uniform mixture of water and surfactant, the system organizes itself in small domains each of which 
has a clear lamellar structure consisting of alternating water-rich and surfactant-rich layers. These domains grow in 
time but not without limit. They are highly dynamic objects which continuously form and disintegrate but whose 
average size does not grow in time once they are formed. In Fig. 13 we display the variation of water and surfactant 
densities within one of the domains, as an example. The densities are averaged over the y— direction within the domain 
and are displayed along the x— axis (the direction of density modulations within this domain). Also shown are d x and 
d y components of surfactant directors. It can be seen that the domain is built up of a stack of surfactant-rich bilayers 
separated by water-rich layers each ~ 2 lattice units wide. The surfactant directors are ordered anti-ferromagneticaly 
within each domain, such that only the hydrophilic heads of surfactant molecules are exposed to water-rich regions. 
In the case of the system with low surfactant density visualization of the data indicates the existence of weak density 
modulations in the system but without any clear domain formation. 

To further quantify the dynamics of self-assembly in the binary water-surfactant system we make use of the 
circularly-averaged surfactant structure factor. In figure 13 we show S(k,t) at timesteps 0, 1000 and 2500 for both 
systems. It can be seen that in both cases S(k,t) has a peak around k — 1.6. This peak corresponds to the repeat 
period of water and surfactant density modulations and its position becomes stable already for t < 100 timesteps. In 
the case of the system with high surfactant concentration we see the emergence of a second peak in S(k,t) at much 
smaller wavectors, indicating that there is another characteristic length-scale in this system, namely the average size 
of lamellar domains. 

In order to identify the driving force behind self-organization of the system we performed two additional calculations, 
using the same initial conditions as before. One simulation was performed with dipolar interactions among surfactant 
molecules switched off (g ss = 0) while in the other simulation we kept g ss = 0.01 but switched off the coupling between 
surfactant and water molecules, by setting g as to 0. In Fig. 14 the spherically averaged structure factor is shown at 
timestep 2500 as obtained from these simulations. For comparison S(k,t) of the full simulations is also shown at the 
same timestep. It can be seen that switching off the interaction between surfactant dipoles result in a structure factor 
which has only a single peak near ft = 1.6 while switching off water-surfactant interaction results in a structure factor 
with only a peak near k = 0. This provide clear evidence that water-surfactant interactions are responsible for the 
formation of alternating water-rich and surfactant-rich layers while formation of domains is a consequence of dipolar 
interactions between surfactant particles. 

IV. CONCLUSIONS 

Building on our recent work JO we developed in this paper a lattice-Boltzmann model for ternary interacting 
amphiphilic fluids. The main features of the model are that interactions among fluid components are realized by 
introducing self-consistently generated mean-field forces between different species and that the orientational degrees 
of freedom of amphiphilic species are explicitly modeled. The mean-field force is incorporated into the scheme in a 
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way which is consistent with the kinetic theory of interacting fluids mixtures and we provided a physical interpretation 
for the action of this force in terms of introducing off-diagonal matrix elements in the BGK collision operator. 

Using a single set of force coupling constants, we have shown that our model exhibits the correct 2D phenomenology 
for both binary and ternary phase systems. Various experimentally observed self-assembling structures form in a 
consistent way as a result of altering the relative amounts of oil, water and amphiphilc in the system. The presence 
of enough surfactant clearly arrests the growth of the domains and we showed that when this happens the final 
domain size is inversely proportional to the amount of surfactant present in the system, in agreement with previous 
molecular dynamics simulations. Our study of the stability of the lamellar structure in 2D confirmed a striking 
difference between lattice-gas and lattice-Boltzmann simulations which result from the absence of fluctuations in 
the lattice-Boltzmann scheme Q]. Self-assembly into local lamellar structure, as found in our simulations, has not 
been reported previously for binary water-surfactant systems but has been observed in microemulsion experiments 
performed on ternary water, oil and surfactant systems p3j] . Our additional calculations indicate that by increasing 
the force coupling constants beyond a certain value the global lamellar phase, which has been reported previously 
JT2| , fl4j , can also be reached in our model. The ability of our model to simulate self-assembly of surfactant aggregates 
in a binary water-surfactant system clearly distinguishes our model from other existing lattice-Boltzmann schemes 
|u] |l6| which do not incorporate the vectorial nature of surfactant molecules and are therefore unable to describe the 
formation of such aggregates. 

Natural refinements of our model would be the inclusion of fluctuations, e.g, via a fluctuating force term compatible 
with the fluctuation-dissipation theorem which enables the system to move out of its metastable states. Also, 
some of the instabilities which we encountered in the present model might be mitigated by using more realistic forms 
for the interaction between different molecules and by adding the Enskog corrections for the collisions in dense systems 
p4j to the BGK collision operator. Our recently developed parallel jl9| version of the lattice-Boltzmann model should 
allow us to extend the present study to 3D for which much more experimental data is available. This will help us to 
select model parameters so that the model will provide a more realistic description of experimental observations; we 
also hope to then study various important phenomena such as flow of complex fluids in porous media. 
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FIG. 1. Snapshots of phase separation of a binary oil-water mixture during a critical quench. From left to right and top 
to bottom timesteps 0, 800, 1600, 2400, 3200, 4800, 6400, 8000, 12000, 16000, 20000, 24000, 30000, , 45000, 60000, 90000 are shown. 
The system size is 256 2 . 



FIG. 2. Temporal evolution of S(k, t) for a critical quench in a binary oil-water system. Timesteps shown are, from bottom 
to top, t = 3200, 4000, 6400, 9600, 12000, 14400, 20000, 24000. The system size is 512 2 . 

FIG. 3. Logarithm of the average domain size R(t) (lattice units) as a function of the logarithm of the time (timesteps) with 
data taken from lattice-Boltzmann simulations of a critical quench in binary immiscible phase separation. The straight line 
corresponds to a growth exponent a = 0.66 and is provided as a guide to the eye only. The system size is 512 2 . 

FIG. 4. Snapshots of time evolution of oil-in-water microemulsion phase. From left to right and top to bottom timesteps 
0, 200, 400, 600, 800, 1200, 1600, 2000, 3000, 4000, 5000, 6000 are shown. The system size is 256 2 . 



FIG. 5. Temporal evolution of the circularly-averaged structure factor S(k, t) for the microemulsion droplet phase shown in 
Figure 4. Timesteps shown, from bottom to top, are t = 200, 1000, 2000, 3000, 5000, 6000. 



FIG. 6. Snapshots of time evolution of sponge microemulsion phase. From left to right and top to bottom timesteps 
0, 200, 400, 600, 800, 1200, 1600, 2000, 3000, 4000, 5000, 6000 are shown. The system size is 256 2 . 



FIG. 7. Temporal evolution of the circularly-averaged structure factor S(k, t) for the sponge microemulsion phase shown in 
Figure 6. Timesteps shown, from bottom to top, are t = 200, 1000, 2000, 3000, 4000, 6000. 

FIG. 8. Time evolution of average domain size R(t) (lattice units) in a ternary system with equal concentrations of water 
and oil (0.5). Curves from top to bottom correspond to systems with average surfactant concentration 0.05,0.1,0.15,0.2,0.3. 
The full lines are the stretched exponential fits to our data. 

FIG. 9. Final domain sizes Rf (lattice units) as a function of the inverse of the surfactant concentration, l/n s . The solid 
line is a linear fit to the first 4 points of our data. 

FIG. 10. Snapshots of the evolution of the lamellar structure in the absence of surfactant. Timesteps shown are clockwise 
from top to bottom t = 0, 300, 600, 3000, 6000, 9000. The system size is 256 2 . 

FIG. 11. Upper panel shows the final state color order parameter, averaged over the y— direction (vertical in Fig. 10) of 
lamellar structure, both with and without surfactant present. The timestep shown is 9000. Note how the presence of surfactant 
sharpens the interfaces between water and oil. Lower panel shows the final state distribution of the surfactant directors, 
averaged over the y— direction at timestep 9000. The system size is 256 2 . 

FIG. 12. Snapshots of self-assembly of a uniform mixture of water and surfactant into lamellar domains. From left to right 
and top to bottom timesteps 0, 300, 3000, 9000. The surfactant-water density difference is shown in grey scaling with black 
corresponding to high surfactant density and white corresponding to high water density. The average concentrations of water 
and surfactant are 0.5 and 0.25, respectively. The system size is 256 2 . 

FIG. 13. Upper panel shows variation of the water density within an "anti-ferromagnetic" domain averaged over the 
y— direction (which is vertical in figure 12) within the domain. Middle panel shows variation of surfactant density within 
the same domain. Lower panel shows the final state distribution of the d x (circles) and d y (squares) of surfactant directors, 
averaged over the y— direction. The timestep shown in all panels is 3000. 
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FIG. 14. Circularly-averaged surfactant density structure factor S(k, t) for binary water and surfactant mixtures shown at 
timesteps and 2500 for a 256 2 binary water-surfactant system. The average density of water is kept at 0.5 while the surfactant 
average density is increased from 0.071 s:w=l:7), to 0.25 (s:w=l:2). For both systems S(k,t) has a peak near k — 1.6 whose 
position corresponds to the repeat period of the lamella. The peak near k = is present only for the system with high surfactant 
concentration and it signals the formation of lamellar domains whose average size correspond to the position of this peak. 

FIG. 15. Circularly-averaged surfactant density structure factor S(k, t) for a binary water and surfactant mixture at timestep 
2500 calculated with both water-surfactant and surfactant-surfactant coupling swtiched on (thick solid line), only surfac- 
tant-water coupling switched on (dashed line) and only surfactant-surfactant coupling switched on (solid line) . The system size 
is 256 2 . The average concentrations of water and surfactant are 0.5 and 0.25 respectively. 
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